Halo_CNV_figure_code_FINAL

library(DNAcopy)
library(foreach)
library(plyr)
library(dplyr)
source("R/CNV_functions_AKS12202017.R")

Introduction

The following code generates figures for the Halobacterium CNV manuscript. This code uses the CNV package DNAcopy to analyze both gene expression and ChIP microarrays.

## All gene expression microarrays (numbers are standardized z scores)
all_data <- read.delim(file = "data/Gene_expression_Master.txt", row.names = 1)
## Empty array to be filled for the generation of segment maps
rarray <- read.delim(file = "data/rarray_allloci.txt")

GE microarray data analysis

The following commented chunk takes a long time to run. A preprocessed object, provided in the github data/ folder, can be loaded instead (see below).

# CNA.ob <- CNA(all_data[,1:1154], all_data$Chr, all_data$Start,data.type = 'logratio', sampleid = colnames(all_data)[1:1154]) 
# smoothed <- smooth.CNA(CNA.ob)
# segsmot <- segment(smoothed, verbose = 0)  ## Set verbose = 1 if number of arrays is not huge
# save(segsmot, file = "data/segmented_GE_arrays.Rdata")

load(file = "data/segmented_GE_arrays.Rdata")

Figure 2, frequency map plots for megaplasmids pNRC100 and pNRC200

Generate maps of segments at different fragment size thresholds

t1f5k.copy <- copymap4(ref_array = rarray, seg_array = segsmot$output, threshold = 1, f.size = 5000)
t1f10k.copy <- copymap4(ref_array = rarray, seg_array = segsmot$output, threshold = 1, f.size = 10000)
t1f20k.copy <- copymap4(ref_array = rarray, seg_array = segsmot$output, threshold = 1, f.size = 20000)
t1f50k.copy <- copymap4(ref_array = rarray, seg_array = segsmot$output, threshold = 1, f.size = 50000)
t1f100k.copy <- copymap4(ref_array = rarray, seg_array = segsmot$output, threshold = 1, f.size = 100000)

Split copy maps by genomic elements

t1f5k.split <- split2(t1f5k.copy)
t1f10k.split <- split2(t1f10k.copy)
t1f20k.split <- split2(t1f20k.copy)
t1f50k.split <- split2(t1f50k.copy)
t1f100k.split <- split2(t1f100k.copy)

pNRC100 by gene expression segment size frequency, figure 2A left

pNRC200 by gene expression segment size frequency, figure 2A right

pNRC100 up and downregulated segment frequency map, Figure 2B left

pNRC200 up and downregulated segment frequency map, Figure 2B right

Figure 3 - GE segmentation up vs downregulation on main chromosome

chromosome, 20kb cutoff, Fig 3A

Calcluate functional enrichment in large GE frequency peaks on chromosome, fig 3B

#Find segments present in >= 5% of arrays, determine functional enrichments in chromosome
ge20k.5<-subset(t1f20k.split[[1]], (Thresh_up + Threshdown) >= 58)
#analysis in this code chunk relies upon https://github.com/amyschmid/histone_arCOG. Load source from there.
GE.5percent.arcog<-cogtest2(ge20k.5$GeneName, cogs, 0.05)
write.table (GE.5percent.arcog, file = "data/ge20k_freqpks_5percent_arCOGs.txt", sep = "\t") #Supplementary table?
#Determine which genes are in which arCOG categories, save as supplementary tables 
ge20k.5.pk1<-cogset(ge20k.5$GeneName, cogs, "Cell motility ")
write.table (ge20k.5.pk1, file = "ge20k_freqpks_5percent_pk1genes.txt",sep = "\t")
ge20k.5.pk2<-cogset(ge20k.5$GeneName, cogs, "Translation; ribosomal structure and biogenesis ")
write.table (ge20k.5.pk2, file = "ge20k_freqpks_5percent_pk2.txt", sep = "\t")
ge20k.5.pk3<-cogset(ge20k.5$GeneName, cogs, "Coenzyme transport and metabolism ")
write.table (ge20k.5.pk3, file = "ge20k_freqpks_5percent_pk3.txt", sep = "\t")
#similar code was used to calculate enrichments in gene expression data frequency peaks for megaplsmids from t1f20k.split[[2]] and t1f20k.split[[3]] except at 5k threshold.

ChIP microarray data analysis

Samples were parsed, median-centered, scaled, and segmented in independent batches before combining. The following code takes segmented chip array data and creates composite maps.

## Read in segmented ChIP data 
test <- read.delim(file = "data/all_clone_fragments_chIP.txt")
## Read in Trmb file 12 for example segmentation output
load(file = "data/trmb12.Rdata")

Example segmentation output, Trmb12 array

shown in Figure 1 example of genomic break point detection in ChIP arrays using DNAcopy

Analyzing: Trmb.12 

Figure 4

Segment grequency map of ChIP CNVs across 48 arrays, categorized by depletion and amplification events

# Generate copymap
crarray <- read.delim(file = "data/rarray_2400_moved_probes.txt")
crarray <- droplevels(crarray[crarray$Chr == "Chr",])
#5k threshold
chip.copy <- copymap4(ref_array = crarray, seg_array = test, threshold = 0.5, f.size = 5000)
chip.split <- split2(chip.copy)[[1]] %>% dplyr::arrange(Start)

Figure 4A

Frequency map of all ChIP CNVs, 5k threshold, 10% array threshold, 16 pks marked.

Fig 4B

Table 1: How many CNVs were detected in each array (GE and ChIP data)?

## First we consider GE data. Filter segments by significance across main chromosome.
all.clone.segs.GE<- segsmot$output
all.clone.segs.GE$seg.length <- (all.clone.segs.GE$loc.end - all.clone.segs.GE$loc.start)
segs.GE.20k <- subset (all.clone.segs.GE, seg.length >=20000)
segs.GE.20k <- subset (segs.GE.20k, chrom == "Chr" & (seg.mean >= 1 | seg.mean <= -1))
#for pNRC200
#segs.GE.20k.p200 <- subset (all.clone.segs.GE, seg.length >=20000)
segs.GE.20k.p200 <- subset (segs.GE.20k.p200, chrom == "pNRC200" & (seg.mean >= 1 | seg.mean <= -1))
#for pNRC100
#segs.GE.20k.p100 <- subset (all.clone.segs.GE, seg.length >=20000)
segs.GE.20k.p100 <- subset (segs.GE.20k.p100, chrom == "pNRC100" & (seg.mean >= 1 | seg.mean <= -1))
## Determine how many significant segments were contained within which GE arrays on main chromosome
freq.per.array20 <- plyr:::count(segs.GE.20k, "ID")
hist(freq.per.array20$freq, col = "grey", main = "Distribution of putative CNVs within each GE array", xlab = "Number of arrays") 
#for pNRC100
freq.per.array20.p100 <- plyr:::count(segs.GE.20k.p100, "ID")
hist(freq.per.array20.p100$freq, col = "grey", main = "Distribution of putative CNVs within each GE array", xlab = "Number of arrays") 

#for pNRC200
freq.per.array20.p200 <- plyr:::count(segs.GE.20k.p200, "ID")
hist(freq.per.array20.p200$freq, col = "grey", main = "Distribution of putative CNVs within each GE array", xlab = "Number of arrays") 

## Next we consider chip data. Filter segments by significance at 5k threshold
test$seg.length <- (test$loc.end - test$loc.start)
segs.chip.5k <- subset (test, seg.length >= 5000)
segs.chip.5k <- subset (segs.chip.5k, chrom == "Chr" & (seg.mean >= 0.5 | seg.mean <= -0.5))
## Determine how many significant segments were contained within which ChIP arrays 
freq.per.array.chip.5 <- plyr:::count(segs.chip.5k, "ID")
hist(freq.per.array.chip.5$freq, col = "grey", main = "Distribution of putative CNVs within each ChIP array", xlab = "Number of arrays")

Figure 5

Frequency map of all ChIP CNVs, with IS elements marked

sig.per.array<-read.delim ("data/num_signif_segs_per_arrays.txt")
hist(sig.per.array$total.signif.segs.per.array, breaks = 12, col = "grey", xlab = "Number CNVs per array", main = "CNVs are frequent")

Fig 5A: IS elements are detected within CNV regions.

Fig 5B: IS elements are associated with CNV regions at higher rate than expected by chance: bootstrapping results

plotchr_chip2(copymap = chip.split, narrays = 48, ymax = 80, ish_table = ish)
legend('topleft',legend=c('CNV','IS', '10% threshold'),col=c('black',rgb(0.1,0.1,.8,.3), 'red'),cex=1.2,bg='white', lty=c(1,2),lwd=2)

abline(h = 10, col = "red")

#points(x = sample(peaks7, 1000), y = c(rep(2,1000)), col = "green", pch =16, cex= .5)

Appendix

hist(ish.boot, col = "gray", xlim = c(0,15), xlab = "Number of IS elements\n in CNV regions", cex.lab=1.2, breaks =10)
abline( v = 14, lty = 2)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.